clear all;

ndraws = 10;
alpha = 0.7; % position-to-position decline in CTR
refctr = [1];
for i=1:11
    refctr = [alpha^i; refctr];
end

randn('seed',3); 

load estrp.mat

postgen = [zeros(size(estrp,1),4)];

for kw = 1:size(estrp,1)
    
    nbids = estrp(kw,1);
    m = estrp(kw,2);
    stdbid = estrp(kw,3);
        
    v = stdbid^2;
    mu = log((m^2)/sqrt(v+m^2));
    sigma = sqrt(log(v/(m^2)+1));
    
    % Generate the profile of random draws
    draws = normrnd(0,1,ndraws,nbids); % generate ndraws random draws
    vals = exp(mu + sigma*sort(draws(:,1:nbids),2));
    
    % The reserve price is either zero if the kw is in the control group or
    % (0.10 + the estimated optimal reserve price)/2 if it's in treatment

    if estrp(kw,4) == 0
        postRP = 0.1;
    else
        postRP = (0.1 + estrp(kw,12))/2;
    end

    % Count the number of bidders above the reserve price for each observation
    bstats = (sum(vals>postRP,2));
                
    % Now iteratively compute equilibrium bids, using the formula from EOS (2007) modified for the case with reserve prices
    vals = max(vals, postRP);
    bids = alpha*postRP + (1-alpha)*vals(:,1);
    for i = 2:(nbids-1)
           bids = [bids alpha*bids(:,i-1)+(1-alpha)*vals(:,i)];
    end
                
    % Now compute the summary statistics for these draws, keeping in mind that
    % observed numbers of bidders will be different for different draws
    mstats =[];
    stdstats =[];
    rps = [];
    for i = 1:ndraws
        if bstats(i) > 1
            mstats = [mstats; mean(bids(i,(nbids-bstats(i)+1):nbids-1))];
            stdstats = [stdstats; std(bids(i,(nbids-bstats(i)+1):nbids-1))];
            if bstats(i) < nbids
                rps = [rps; estrp(kw,7) * (bids(i,(nbids-bstats(i)):nbids-1) * refctr(13-bstats(i):12))];
            else
                rps = [rps; estrp(kw,7) * ([postRP bids(i,(nbids-bstats(i)+1):nbids-1)] * refctr(13-bstats(i):12))];
            end
        elseif bstats(i) == 1
            rps = [rps; postRP*estrp(kw,7)];
        elseif bstats(i) == 0
            rps = [rps; 0];
        end
    end
    
    postgen(kw,:) = [mean(bstats) mean(mstats) mean(stdstats) mean(rps)];

end

finalDataSet = [estrp postgen];

% Save the synthetic dataset
writematrix(finalDataSet);
save finalDataSet.mat finalDataSet
